%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plots for simulated data
% this version: October 2013
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear; 

% load data for single plots
prodfnc = 5;

if prodfnc == 1
    load('response_067_not')
    load('theorirf_phi_min495')
    theorirf = theorirf_min495;
elseif prodfnc == 2
    load('response_phi_145')     
    load('theorirf_phi_145')
    theorirf = theorirf_145;
elseif prodfnc == 3
    load('response_267_not')
    load('theorirf_phi_539')
    theorirf = theorirf_539;
elseif prodfnc == 4
    load('response_phi_625')
    load('theorirf_phi_625')
    theorirf = theorirf_625;
elseif prodfnc == 5
    load('response_167_not')
    load('theorirf_phi_401')
    theorirf = theorirf_401;
elseif prodfnc == 6
    load('response_phi_685')
    load('theorirf_phi_685')
    theorirf = theorirf_685;
elseif prodfnc == 7
    load('response_500_not')
    load('theorirf_phi_800')
    theorirf = theorirf_800;
elseif prodfnc == 8
    load('response_phi_900')
    load('theorirf_phi_900')
    theorirf = theorirf_900;
elseif prodfnc == 9
    load('response_phi_800_dc7')
    load('theorirf_phi_800')
elseif prodfnc == 10
    load('response_phi_900_dc7')
    load('theorirf_phi_900')
elseif prodfnc == 11
    load('response_phi_min_495_dc7')
    load('theorirf_phi_min495')
    theorirf = theorirf_min495;
end

Respmedall = Respmedall*100;
simsize = size(Respmedall,3);
for h=1:size(Respmedall,2)
    Respmed(:,h) = sum(Respmedall(:,h,:),3)/simsize;
    for n=1:size(Respmedall,1)
        Respstd(n,h) = std(Respmedall(n,h,:))/sqrt(simsize);
        Resplow(n,h) = Respmed(n,h)-Respstd(n,h);
        Resphigh(n,h) = Respmed(n,h)+Respstd(n,h);
    end
end

% normalization
MAXits=2000;
tol=1.0e-6;
% to match 20qtr values (from Fisherresponse)
iresp = 0.9361;%1.0076;
nresp = 0.3898;%0.3513;
iscalar_Respmed = Respmed(1,10)/iresp;
nscalar_Respmed = Respmed(6,10)/nresp;
iscalar_theory = theorirf(10,1,2)/(-iresp);
nscalar_theory = theorirf(10,2,1)/nresp;

do_normalize = 1;
if do_normalize == 1
    Respmed(1:4,:) = Respmed(1:4,:)/iscalar_Respmed;
    Resplow(1:4,:) = Resplow(1:4,:)/iscalar_Respmed;
    Resphigh(1:4,:) = Resphigh(1:4,:)/iscalar_Respmed;
    Respmed(5:8,:) = Respmed(5:8,:)/nscalar_Respmed;
    Resplow(5:8,:) = Resplow(5:8,:)/nscalar_Respmed;
    Resphigh(5:8,:) = Resphigh(5:8,:)/nscalar_Respmed;
    
    theorirf(:,:,2) = theorirf(:,:,2)/iscalar_theory;
    theorirf(:,:,1) = theorirf(:,:,1)/nscalar_theory;
end

% saving the data
if do_normalize == 1
    if prodfnc == 1
        Respmed_min495 = Respmed;
        theorirf_min495 = theorirf;
        save Respmed_phi_min495 Respmed_min495
        save theorirf_phi_min495 theorirf_min495
    elseif prodfnc == 2
        Respmed_145 = Respmed;
        theorirf_145 = theorirf;
        save Respmed_phi_145 Respmed_145
        save theorirf_phi_145 theorirf_145
    elseif prodfnc == 3
        Respmed_539 = Respmed;
        theorirf_539 = theorirf;
        save Respmed_phi_539 Respmed_539
        save theorirf_phi_539 theorirf_539
    elseif prodfnc == 4
        Respmed_625 = Respmed;
        theorirf_625 = theorirf;
        save Respmed_phi_625 Respmed_625
        save theorirf_phi_625 theorirf_625
    elseif prodfnc == 5
        Respmed_401 = Respmed;
        theorirf_401 = theorirf;
        save Respmed_phi_401 Respmed_401
        save theorirf_phi_401 theorirf_401
    elseif prodfnc == 6
        Respmed_685 = Respmed;
        theorirf_685 = theorirf;
        save Respmed_phi_685 Respmed_685
        save theorirf_phi_685 theorirf_685
    elseif prodfnc == 7
        Respmed_800 = Respmed;
        theorirf_800 = theorirf;
        save Respmed_phi_800 Respmed_800
        save theorirf_phi_800 theorirf_800
    elseif prodfnc == 8
        Respmed_900 = Respmed;
        theorirf_900 = theorirf;
        save Respmed_phi_900 Respmed_900
        save theorirf_phi_900 theorirf_900
    elseif prodfnc == 9
        Respmed_800_dc7 = Respmed;
        save Respmed_phi_800_dc7 Respmed_800_dc7
    elseif prodfnc == 10
        Respmed_900_dc7 = Respmed;
        save Respmed_phi_900_dc7 Respmed_900_dc7
    elseif prodfnc == 11
        Respmed_min_495_dc7 = Respmed;
        save Respmed_phi_min_495_dc7 Respmed_min_495_dc7
    end
else
    if prodfnc == 1
        Respmed_min495 = Respmed;
        save Respmed_phi_min495_nn Respmed_min495
    elseif prodfnc == 2
        Respmed_145 = Respmed;
        save Respmed_phi_145_nn Respmed_145
    elseif prodfnc == 3
        Respmed_539 = Respmed;
        save Respmed_phi_539_nn Respmed_539
    elseif prodfnc == 4
        Respmed_625 = Respmed;
        save Respmed_phi_625_nn Respmed_625
    elseif prodfnc == 5
        Respmed_401 = Respmed;
        save Respmed_phi_401_nn Respmed_401
    elseif prodfnc == 6
        Respmed_685 = Respmed;
        save Respmed_phi_685_nn Respmed_685
    elseif prodfnc == 7
        Respmed_800 = Respmed;
        save Respmed_phi_800_nn Respmed_800
    elseif prodfnc == 8
        Respmed_900 = Respmed;
        save Respmed_phi_900_nn Respmed_900
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plots

HORIZON_impresp = 10;

Respmed = Respmed(:,1:HORIZON_impresp);
Resplow = Resplow(:,1:HORIZON_impresp);
Resphigh = Resphigh(:,1:HORIZON_impresp);
        
figure(1);
tt = (0:HORIZON_impresp);
zerobar = zeros(1,HORIZON_impresp);
subplot(2,4,1); plot(tt,[0,zerobar],'-k',tt,[nan,-Respmed(1,:)],'-k',tt,[nan,-Resplow(1,:)],'--r',tt,[nan,-Resphigh(1,:)],'--r',tt,[nan,theorirf(2:HORIZON_impresp+1,1,2)'],'--k')
title('price');
ylabel('IBT shock');
if prodfnc == 5
    axis([0 10 -1.5 0.1]);
end
subplot(2,4,2); plot(tt,[0,zerobar],'-k',tt,[nan,Respmed(2,:)],'-k',tt,[nan,Resplow(2,:)],'--r',tt,[nan,Resphigh(2,:)],'--r',tt,[nan,theorirf(2:HORIZON_impresp+1,2,2)'],'--k')
title('prod.')
if prodfnc == 5
    axis([0 10 -0.1 0.5]);
end
subplot(2,4,3); plot(tt,[0,zerobar],'-k',tt,[nan,Respmed(4,:)],'-k',tt,[nan,Resplow(4,:)],'--r',tt,[nan,Resphigh(4,:)],'--r',tt,[nan,theorirf(2:HORIZON_impresp+1,3,2)'],'--k')
title('tot. hours');
if prodfnc == 5
    axis([0 10 -0.1 0.5]);
end
subplot(2,4,4); plot(tt,[0,zerobar],'-k',tt,[nan,Respmed(3,:)],'-k',tt,[nan,Resplow(3,:)],'--r',tt,[nan,Resphigh(3,:)],'--r',tt,[nan,theorirf(2:HORIZON_impresp+1,4,2)'],'--k')
title('premium');
if prodfnc == 5
    axis([0 10 -0.1 0.1]);
end
subplot(2,4,5); plot(tt,[0,zerobar],'-k',tt,[nan,-Respmed(5,:)],'-k',tt,[nan,-Resplow(5,:)],'--r',tt,[nan,-Resphigh(5,:)],'--r',tt,[nan,theorirf(2:HORIZON_impresp+1,1,1)'],'--k')
ylabel('INT shock');
xlabel('Quarters');
if prodfnc == 5
    axis([0 10 -1.5 0.1]);
end
subplot(2,4,6); plot(tt,[0,zerobar],'-k',tt,[nan,Respmed(6,:)],'-k',tt,[nan,Resplow(6,:)],'--r',tt,[nan,Resphigh(6,:)],'--r',tt,[nan,theorirf(2:HORIZON_impresp+1,2,1)'],'--k')
xlabel('Quarters');
if prodfnc == 5
    axis([0 10 -0.1 0.5]);
end
subplot(2,4,7); plot(tt,[0,zerobar],'-k',tt,[nan,Respmed(8,:)],'-k',tt,[nan,Resplow(8,:)],'--r',tt,[nan,Resphigh(8,:)],'--r',tt,[nan,theorirf(2:HORIZON_impresp+1,3,1)'],'--k')
xlabel('Quarters');
if prodfnc == 5
    axis([0 10 -0.1 0.5]);
end
subplot(2,4,8); plot(tt,[0,zerobar],'-k',tt,[nan,Respmed(7,:)],'-k',tt,[nan,Resplow(7,:)],'--r',tt,[nan,Resphigh(7,:)],'--r',tt,[nan,theorirf(2:HORIZON_impresp+1,4,1)'],'--k')
xlabel('Quarters');
if prodfnc == 5
    axis([0 10 -0.1 0.1]);
end
if prodfnc == 1
    set(gcf,'paperunits','centimeters','PaperPosition',[0 0 13 8])
    print -depsc2 model_simulirfs_1_0511.eps
elseif prodfnc == 3
    set(gcf,'paperunits','centimeters','PaperPosition',[0 0 13 8])
    print -depsc2 model_simulirfs_3_0511.eps
elseif prodfnc == 5
    set(gcf,'paperunits','centimeters','PaperPosition',[0 0 13 8])
    print -depsc2 model_simulirfs_5_0511.eps
elseif prodfnc == 7
    set(gcf,'paperunits','centimeters','PaperPosition',[0 0 13 8])
    print -depsc2 model_simulirfs_5_0511.eps
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
do_normalize = 1;
% joint plot
load('Response_spec5_sample2_prem_1_prod0_price0_hours3_break00_est183.mat')
% sorting the Response
for j=1:(size(Response,1))
    Respsort(j,:,:) = sort(Response(j,:,:),3);
end;

HORIZON_ident = 80;
ndraws = size(Response,3);
% Bayesian one-standard error
FRespmed(:,:) = Respsort(:,1:HORIZON_ident,fix(0.5*ndraws));
FResplow(:,:) = Respsort(:,1:HORIZON_ident,fix(0.16*ndraws));
FResphigh(:,:) = Respsort(:,1:HORIZON_ident,fix(0.84*ndraws));

if do_normalize == 1
    load('theorirf_phi_min495')
    load('theorirf_phi_145')
    load('theorirf_phi_539')
    load('theorirf_phi_625')
    load('theorirf_phi_800')
    load('theorirf_phi_900')

    HORIZON_impresp = 10;
    figure(2)
    tt = (0:HORIZON_impresp);
    zerobar = zeros(1,HORIZON_impresp);
    plot(tt,[0,zerobar],'-k',tt,[nan,FRespmed(3,1:HORIZON_impresp)],'-k',tt,[nan,FResplow(3,1:HORIZON_impresp)],':r',tt,[nan,FResphigh(3,1:HORIZON_impresp)],':r',tt,[nan,theorirf_min495(2:HORIZON_impresp+1,4,2)'],'--k',tt,[nan,theorirf_145(2:HORIZON_impresp+1,4,2)'],'--k',tt,[nan,theorirf_539(2:HORIZON_impresp+1,4,2)'],'--k',tt,[nan,theorirf_625(2:HORIZON_impresp+1,4,2)'],'--k',tt,[nan,theorirf_800(2:HORIZON_impresp+1,4,2)'],'--k')
    title('premium');
    ylabel('IBT shock');
    xlabel('Quarters');
    set(gcf,'paperunits','centimeters','PaperPosition',[0 0 10 6])
    print -depsc2 model_comparison_1211.eps


else
    load('Respmed_phi_min495_nn')
    load('Respmed_phi_145_nn')
    load('Respmed_phi_539_nn')
    load('Respmed_phi_625_nn')
    load('Respmed_phi_800_nn')
    load('Respmed_phi_900_nn')

    HORIZON_impresp = 50;
    figure(2)
    tt = (0:HORIZON_impresp);
    zerobar = zeros(1,HORIZON_impresp);
    plot(tt,[0,zerobar],'-k',tt,[nan,FRespmed(9,1:HORIZON_impresp)],'-k',tt,[nan,FResplow(9,1:HORIZON_impresp)],'--r',tt,[nan,FResphigh(9,1:HORIZON_impresp)],'--r',tt,[nan,Respmed_min495(3,1:HORIZON_impresp)],'-k',tt,[nan,Respmed_145(3,1:HORIZON_impresp)],'-k',tt,[nan,Respmed_539(3,1:HORIZON_impresp)],'-k',tt,[nan,Respmed_625(3,1:HORIZON_impresp)],'-k',tt,[nan,Respmed_800(3,1:HORIZON_impresp)],'-k')
    title('premium');
    ylabel('invest. shock');
    xlabel('Quarters');

end